/*
 <one line to give the library's name and an idea of what it does.>
 Copyright (C) 2013  Guillaume <email>
 
 This library is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation; either
 version 2.1 of the License, or (at your option) any later version.
 
 This library is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public
 License along with this library; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */


#ifndef MSKINETICOPERATOR_H
#define MSKINETICOPERATOR_H

#include <msScalarFunction.h>
#include <msSystem.h>

namespace atomism
{
    
    
    class msLinearScalarFunction : public msScalarFunction
    {
        
        static msRegistrar Registrar;
        static bool  isLinearScalarFunctionRegisteredInPython;
        static msRegister* createInstance() { return new msLinearScalarFunction; }
        
    protected:
        
        void registryInPython();
        
        void initialize()  {
            msScalarFunction::initialize();
            msTreeMapper::declareAttribute(a,"a");
            msTreeMapper::declareAttribute(b,"b");
            msTreeMapper::declareAttribute(Indice,"Indice");
        };
        
        void update(){  msScalarFunction::update(); }
        
    public:
        
        msLinearScalarFunction() : msScalarFunction(){
            
            constructVar("msLinearScalarFunction","msLinearScalarFunction","linear scalar function");
        }
        
        static boost::shared_ptr<msLinearScalarFunction> New(string unitexpr ){
            
            boost::shared_ptr<msLinearScalarFunction> T( new msLinearScalarFunction );
            LOGGER_ENTER_FUNCTION_DBG("msLinearScalarFunction::New(boost::shared_ptr<msUnitsManager> unit)",T->getFullId());
            T->initialize();
            T->getYUnit().set(unitexpr);
            T->update();
            LOGGER_EXIT_FUNCTION();
            return T;
        }
        
        virtual boost::shared_ptr<msTreeMapper> clone() { boost::shared_ptr<msLinearScalarFunction> clone = New("");
            clone->msTreeMapper::operator=(*this);
            return clone;
        }
        
        boost::shared_ptr<msTreeMapper> setCoefsIndice(double a0, double a1,int i);
        
        double evaluate();
        
      	std::string getLabel(){ return ("Linear scalar function"); }
        
        ostream&  print(ostream& out) const {  msScalarFunction::print(out);
            output::printHeader(out,"Linear scalar function")<<"\n a="<<a<<", b="<<b<<endl;
            return out;
        }
        
    private:
        
        double a; double b;
        int Indice;
    };
    
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------

    class msDynamicDof : public msTreeMapper
    {
        friend class msKineticOperator;
        
    private:
        //! @Name from msRegistrer
        //@{
        static msRegistrar Registrar;
        static msRegister* createInstance() { return new msDynamicDof; }
        static bool isDynamicDofRegisteredInPython;
        //@}
        
    protected:
        
        void initialize(){msTreeMapper::initialize();
            declareChild<msScalarFunction>( Function , msScalarFunction::New(""), "Function" );
            declareChild<msEntity>( Entity , msEntity::New(msUnitsManager::New("")),   "Entity" );
            declareAttribute( Indice, "Indice");
        }
        
        void update(){
            msTreeMapper::update();
        }
        msDynamicDof():msTreeMapper() {
            constructVar("msDynamicDof","DynamicDof","dynamic variable");
        }
        
        void registryInPython();
        
    public:
        
        static boost::shared_ptr<msDynamicDof> New(){
            
            boost::shared_ptr<msDynamicDof> T( new msDynamicDof );
            LOGGER_ENTER_FUNCTION_DBG("msDynamicDof::New(boost::shared_ptr<msUnitsManager> units)",T->getFullId());
            T->initialize(); T->update();
            LOGGER_EXIT_FUNCTION();
            return T;
        }
        
        virtual boost::shared_ptr<msTreeMapper> clone() {
            boost::shared_ptr<msDynamicDof> clone = New();
            clone->msTreeMapper::operator=(*this);
            return clone;
        }
        
        void set(boost::shared_ptr<msScalarFunction> function,
                 boost::shared_ptr<msEntity> entity, size_t i);
        
        std::ostream& print(std::ostream& out) const {
            
            out<<" Degree of freedom no. "<<Indice<<" for the entity of id '"<<Entity->getId()
            <<"'; scalar function:"<<Function->getType();
            return out;
        }
        
    private:
        
        msChild<msScalarFunction>   Function;
        msChild<msEntity>           Entity;
        size_t Indice;
    };
    
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
   
    class msKineticOperator;
    
    class msJacobianOfDispl
    {
        friend class msSystem;
        friend class msKineticOperator;
        
    public:
        
        msJacobianOfDispl();
        
        void initialize( msKineticOperator* kinetic );
        
        void compute( const msUnitsManager& units, bool isolated );
        
        void resize(size_t noOfat,size_t noOfdof);
        
        const Eigen::Vector3d& operator()(int elem,int qi);
        
        void setJacColumn( const msGeneralizedCoordinates::UnfreezIt q , size_t ind, const msUnitsManager& units );
        
        void setDisplacment( const msVector3dContainer& coors0,
                            const msVector3dContainer& coors1,
                            size_t ind, double conv );
        
        std::ostream& print(std::ostream& out) const;
        
    private:
        
        boost::numeric::ublas::matrix<Eigen::Vector3d> Displ_Elem_qi;
        
        msKineticOperator* 	      KineticOperator;
        msSystem* 	 	      System;
        Eigen::Vector3d*     	      OverallTrans;
        Eigen::Matrix3d*     	      OverallRot;
        msGeneralizedCoordinates*       GeneralizedCoordinates;
        
        size_t NoOfAtoms, NoOfDofs;
        
        msVector3dContainer 	      Coors0;
        Eigen::Vector3d 		      CenterOfMass0;
        msVectorOfDofs 		      Dofs0;
        
        msVector3dContainer 	      Coors1;
        Eigen::Vector3d 		      CenterOfMass1;
        msVectorOfDofs		      dDofs;
    };
    
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    //-----------------------------------------------------------------------------
    
    class msKineticOperator : public msPhysicalInterface
    {
        friend class msJacobianOfDispl;
        
        class msKineticOperatorParams : public msParamsManager
        {
            static msRegistrar Registrar;
            static msRegister* createInstance() { return new msKineticOperatorParams; }
            
        public:
            msKineticOperatorParams() : msParamsManager() {
                /* addParameter( "MaxStep" , "25" ,   "Max. steps for J annihilation" , msUnit() );
                 addParameter( "ResEkin" , "1e-5" , "Max. steps for J annihilation" , msUnit() );*/
            }
        };
    private:
        //! @Name from msRegistrer
        //@{
        static msRegistrar Registrar;
        static bool  isKineticOperatorRegisteredInPython;
        static msRegister* createInstance() { return new msKineticOperator; }
        //@}
        
    protected:
        
        void registryInPython();
        
        void initialize(){msPhysicalInterface::initialize();
			declareChildren<msDynamicDof>( DynamicDofs, "DynamicDofs" );
			declareChild<msSystem>( System , msSystem::New(getUnits()), "System" );
			declareChild<msGeneralizedCoordinates>( GeneralizedCoordinates ,
                                                   msGeneralizedCoordinates::New(getUnits()),
                                                   "GeneralizedCoordinates" );
        }
        void update(){msPhysicalInterface::update();JacobianOfDispl.reset(new msJacobianOfDispl);
            KMatrix.reset(new msKineticMatrix(GeneralizedCoordinates->noOfActive(),
		                                      GeneralizedCoordinates->noOfActive()));
        }
        
        msKineticOperator():msPhysicalInterface() {  constructVar("msKineticOperator","KineticOperator","kinetic operator");
        }
    public:
        
        enum TestType { system , coordinates };
        
        typedef Eigen::MatrixXd msKineticMatrix;
        
        static boost::shared_ptr<msKineticOperator> New( boost::shared_ptr<msUnitsManager> units ){
            
            boost::shared_ptr<msKineticOperator> T( new msKineticOperator );
            LOGGER_ENTER_FUNCTION_DBG("msKineticOperator::New(boost::shared_ptr<msUnitsManager> units)",T->getFullId());
            T->initialize();
            T->setUnits(units);
            T->setParameters( boost::shared_ptr<msParamsManager> (new msKineticOperatorParams()) );
            T->getParameters()->setAffiliation( T->mySharedPtr());
            T->update();
            LOGGER_EXIT_FUNCTION();
            return T;
	    }
        virtual boost::shared_ptr<msTreeMapper> clone() {
            boost::shared_ptr<msKineticOperator> clone = New( msUnitsManager::New() );
            clone->msTreeMapper::operator=(*this);
            return clone;
        }
        /* msKineticOperator(const msKineticOperator& rhs);
         
         msKineticOperator& operator=(const msKineticOperator& rhs);*/
        
        boost::shared_ptr<msTreeMapper> set( boost::shared_ptr<msGeneralizedCoordinates> coords ,
                                            boost::shared_ptr<msSystem> system );
        
        boost::shared_ptr<msTreeMapper> setCoordinates( boost::shared_ptr<msGeneralizedCoordinates> coords);
        boost::shared_ptr<msTreeMapper> setSystem( boost::shared_ptr<msSystem> system );
        
        boost::shared_ptr<msTreeMapper> addDynamicDof( std::string id , size_t indice ,
                                                      boost::shared_ptr<msScalarFunction> fct);
        
        boost::shared_ptr<msTreeMapper> addDynamicDof( boost::shared_ptr<msEntity> entity , size_t indice ,
                                                      boost::shared_ptr<msScalarFunction> fct);
        
        boost::shared_ptr<msTreeMapper> addDynamicDofLinear( boost::shared_ptr<msEntity> entity , size_t indice ,
                                                            boost::shared_ptr<msScalarVariable> q );
        
        void setDynamicDofs();
        
        void setPositions(){System->computeCartCoordinates();};
        
        msKineticWork move(  size_t i );
        
        double kineticFunction(){ return kineticFunction(*getUnits()); };
        
        double kineticFunction(const msUnitsManager& units);
        
        const msKineticMatrix&   computeKMat() { return computeKMat(*getUnits());}
        
        const msKineticMatrix&   computeKMat(const msUnitsManager& units);
        
        //    void computeDofsVariation( msScalarVariable& qi );
        
        const msKineticMatrix& getKMatrix() { return *KMatrix; }
        
        const boost::shared_ptr<msGeneralizedCoordinates> getCoordinates() const {
            return GeneralizedCoordinates.getSharedPtr(); }
        
        boost::shared_ptr<msGeneralizedCoordinates>       getCoordinates() {
            return GeneralizedCoordinates.getSharedPtr(); }
        
        boost::shared_ptr<msSystem> getSystem(){ return System.getSharedPtr();}
        
        std::ostream& print(std::ostream& out) const;
        
    private:
        
        void registerMembers();
        
        msChildren<msDynamicDof> 		   DynamicDofs;
        msChild<msSystem> 		   System;
        msChild<msGeneralizedCoordinates>    GeneralizedCoordinates;
        
        
        boost::shared_ptr<msKineticMatrix>   KMatrix;
        boost::shared_ptr<msJacobianOfDispl> JacobianOfDispl;
        
        msVectorOfDofs dDofs,Dofs0;
        msVector3dContainer Coors0,Coors1;
        
        bool checkTest(TestType test) const;    
    };
    
    ostream& operator<<(ostream&,const msKineticOperator&);
    ostream& operator<<(ostream&,const msJacobianOfDispl&);
}
#endif // MSKINETICOPERATOR_H
